Statistics of Deep Energy States in Coulomb Glasses 
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We study the statistics of local energy minima in the configuration space and the energy relax- 
ation due to activated hopping in a system of interacting electrons in a random environment. The 
distribution of the local minima is exponential, which is in agreement with extreme value statistics 
considerations. The relaxation of the system energy shows logarithmic time dependence reflecting 
the ultrametric structure of the system. 



Possible phases and the corresponding dynamic behav- 
iors of strongly correlated disordered systems remain one 
of the central unresolved issues. Doped semiconductors 
in the insulating state are an exemplary system endowed 
with strong long-range Coulomb interactions and strong 
disorder. It was hypothesized that the combined action 
of both inflict the glassy phase, resulting in the term 
electron- or Coulomb glass Ever since, a Coulomb 
glass has been a subject of the intense research, see Q; 
however, despite the impressive developments, a consis- 
tent picture of its nature is still not available and the 
consensus on the interrelation between its inherent glassy 
ingredients like Coulomb gap, slow dynamics and aging, 
the nature of the low-lying metastable states and the 
mechanisms of hopping transport is not yet achieved. A 
classic qualitative derivation [3[ of the Coulomb gap in- 
duced in the single-particle density of states due to pair 
correlations remains the top result in the field. 

One of the main characteristics of the glassy state is the 
existence of an infinite number of low-lying states (val- 
leys in the rugged free energy relief) separated by barriers 
growing infinitely in the thermodynamic limit. This pic- 
ture emerged from the study of the infinite range interac- 
tion model for spin glasses (Sherrington-Kirkpatrick spin 
glass) Q and was later shown to imply an ultrametric 
structure of the configurational space |5|. 

There was important recent progress in relating the 
appearance of the Coulomb gap in electron systems to 
the generic glassy behavior [6[ , where a lattice model of 
spinless interacting electrons in the presence of random- 
ness was considered and a nonlinear screening theory was 
formulated to tackle Coulomb gap formation and its rela- 
tion to glassy freezing. The formation of the glass phase 
with a large number of metastable states was identified 
via the emergence of a replica symmetry breaking insta- 
bility. In [7[, the connection between the glass phase 
characteristics and a Coulomb gap was argued, notic- 
ing that within the locator-approximation the correlated 
electron system maps to the Sherrington-Kirkpatrick spin 
glass. Remarkably, recent experimental studies of low- 
frequency resistance noise in silicon inversion layers in- 
dicated glassy freezing and the presence of long-range 



correlations, consistent with the hierarchical picture of 
glassy dynamics This gives strong support to a hi- 
erarchical structure of experimentally observed Coulomb 
glasses and calls for further comprehensive study of hier- 
archical glassy systems. 

In this Letter we demonstrate that the distribution of 
low-lying energy states in systems with an ultrametric 
configurational space obeys extreme value statistics giv- 
ing rise to an exponential distribution for the extremely 
low energies, p oc exp(—E/Eo), where Eq is a character- 
istic energy. The exponential tail is verified numerically 
for the lattice model of a disordered Coulomb system. 
We find that system's energy relaxes with time following 
the logt law, as indeed expected for ultrametric hier- 
archical structures This offers strong evidence that 
electronic systems with long-range Coulomb interactions 
form a Coulomb glass with exponential distribution of 
lowest energy states. 

Exponential distribution. — We start with the deriva- 
tion of the exponential statistics for the low-lying energy 
states of systems with ultrametric topology of the con- 
figurational space. This space can be represented graph- 
ically as a hierarchical tree which is the union (or sum) 
of hierarchical trees with branching numbers n G N\{1} 
(n- trees, see Fig. [I]). The crown of the tree (leaves) cor- 
responds to a set of nearly equal energies that form an 
almost degenerate ground state of the system. We con- 
sider now the subset of the configurational space charac- 
terized by an n-tree and find the distribution function of 
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FIG. 1: Illustration of a 3-tree, an element of the hierarchical 
tree representing the configurational space with ultrametric 
topology. 
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the corresponding lowest energy states. It is convenient 
to introduce positive variables {E} measuring the depths 
of the potential wells. Starting with the lowest hierarchi- 
cal level (crown) consider one of the sibling states sets, 
{E\° , i G N- n }, where the superscript (0) marks the low- 
est hierarchy level. We define = maxi<i<„{£' J - ' ) } 
and assign it to the branching point level (1) generat- 
ing the considered 0-th set of siblings. Now we have a 
collection of level-(l) sibling n-sets. Repeating the proce- 
dure, we obtain the next generation of low-lying energies 
{E^ } that constitute the assortment of the second level 
sibling n-sets. The described iteration should converge to 
a limiting distribution for the probability V n (E^) for 
the lowest states at the j-level of the hierarchy, as j — * oo. 
This distribution must be stable under the max opera- 
tion. In other words if E = max; {£".;} then the proba- 
bility distribution of the extrema E is the same as that 
of each member Ei of the set. The distribution we seek 
belongs to the class of so-called extreme statistics [l(J ■ 

To briefly summarize the theory of extreme distribu- 
tions, we consider a set of identically distributed indepen- 
dent random variables Xi , with 1 < i < n (in the present 
problem these are all the low-lying energies within one 
set of sibling states). Let M n = m&xi<i< n {Xi} and let 
V\ (x) be the probability that any of the random variables 
Xi is less than x. The probability that the largest value, 
M n , of the set {Xi;i <E N-' 1 } is less than x is simply the 
probability that all of the XiS fall short of x. Since the 
variables are independent, this is given by 



V [ max {Xi} < x 

KKn 



V n (x) = [P 1 (x)] n . (1) 



The central concept of the theory of extreme statistics 
is the stability postulate (Io| , stating that after many it- 
erations the distribution describing M n for a sufficiently 
large set {Xi} is the same as the distribution for each 
variable Xi. This requires that the probability V n (x) sat- 
isfied the functional equation V n (x) = [P„(a n x + b n )] n , 
where a n and b„ arc positive constants. It can be shown 
that, under rather nonrestrictive assumptions about the 
initial distribution of random variables {Xi} there are 
only three classes of distribution functions depending on 
its support. In particular, if x £ R, 

V n {x) = exp[-ncxp(-x)] , a n = 1, b n = ln(n), (2) 

and, therefore, the limiting form for a distribution of the 
lowest energy states corresponding to an n-tree subset of 
the configurational space is V n (E) = exp[— nexp(— E)], 
where E is measured in dimensionless energy units. 

As mentioned before, the configurational space of our 
system endows with an ultrametric structure which can 
be represented as an union of all possible n-trees. There- 
fore the global distribution density can be found as 



where w(n) is the weight of the n-trees in the config- 
urational space. The weight w(n) is determined by the 
physical properties of the system in question. In our case 
one observes that the elemental process of generating two 
nearly degenerate energy states is the pair electron ex- 
chan ge w ithin a specially arranged compact 4-sites aggre- 
gate [11| . Analogously, compact aggregates allowing the 
exchange of n pairs, generate the same number of sibling 
states of the corresponding hierarchical trees. Thus w(n) 
is nothing but the percolation cluster size-distribution 
function: w n oc n _s , where ,s > 2. Plugging this weight 
into Eq. Q one arrives at the exponential density distri- 
bution for the low-lying energies, p(E) oc exp(— E/Eq), 
where s is absorbed into the normalization energy Eq. 

The extreme value distributions for extremely low- 
energy states belongs to the so-called were derived [12| for 
several systems and problems possessing one step replica 
symmetry breaking solutions, such as the random energy 
model and decaying Burgers turbulence. It was found 
that extremely low energy states of elastic manifolds in 
random environment (modeling a rich variety of physi- 
cal systems including vortices in superconductors, charge 
density waves, Wigner crystals, domain walls, and many 
others) also obey exponential statistics resulting in creep 
motion which is specific for glassy states Our re- 

sults taken together with those of allow to ar- 

rive at the conclusion that the exponential distribution 
of the extremely low-lying energy states, the resulting 
creep response, and the related 1// noise [LJ] are generic 
properties of all strongly correlated disordered systems 
possessing a fractal structure of configurational space. 

Model. — Now we turn to numerical studies of the en- 
ergy configurational space. The lattice Coulomb model 
was studied numerically in the past 15, [Til . 17 1, the 
typical system size was about 10 2 . Optimized numer- 
ical procedures enabled us to consider 100-times larger 
systems and analyze about 10 5 different energy minima. 
The Hamiltonian of a 2D Coulomb lattice with disorder 




p{E) = J dnw(n) 



dV„(E) 
dE ' 



(3) 



FIG. 2: Energy distribution of about 12000 local minimum 
states in half-log representation with linear fit. Linear scale 
with exponential fit (triangular inset). 
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is 

H = Y,^ e i = Ua i n i + -22— — " . ( 4 ) 

where i and j label the lattice sites, = | — r j I > 
U Ui with on G [— 1; 1] define the quenched uniformly dis- 
tributed random site energies, and rii are the occupation 
numbers (values or 1) with average v. 

The simulations were performed on the lattices with 
typically N = 100 2 sites with periodic boundary condi- 
tions and the filling factor v = 0.5. We took the disorder 
spread, U , be equal to the Coulomb interaction at the dis- 
tance of the lattice constant. Consequently, U, the lattice 
constant, and e were set to unity; vectors label square 
lattice sites. The behavior of the system is governed by 
the interplay between disorder and Coulomb interactions. 
At very weak disorder the system is crystal-like; increase 
in U brings glassy behavior. 

Local minima. — To track down local minima of the 
system free energy, we choose an arbitrary initial dis- 
tribution of electrons and perform sequential transfer of 
electrons between the adjacent sites (of course such an 
elemental hop can be realized only if the target neigh- 
boring site is vacant). At each step the total Coulomb 
energy of the system is calculated using Eq.(0|. If any 
nearest site hop gives rise to the increase in total en- 
ergy, as compared to its value at previous step, the state 
and the corresponding energy is marked as energy min- 
imum. In subsequent iterations one finds transfers that 
decrease total energy executing thus dissfusion of the sys- 
tem into more and more deep energy valleys. Note that 
this algorithm lists all the possible low-lying states. To 
ensure evolving the system into true low lying minima 
we were starting collection statistics after reaching the 
states where at least 10 5 random hops were needed to 
go from one local minimum to another. Shown in Fig. [5] 
is the obtained distribution of energies corresponding to 
local minima; one sees that the distribution tail obeys 
the exponential law, N min (E) ~ e~ E/E °. 

For all the computed minima we define the normalized 
site occupation number difference A a p = N^ 1 \nf — 
nf | corresponding to all pairs of minima {a, (3}. The dis- 
tribution of A a p, i. e., the fraction of the system which 
has different occupation numbers for two different min- 
ima is shown in Fig. [3J To drive our system (the filling 
factor 0.5) from one minimum to another, about 5.4% 
of the electrons has to be transferred to different lattice 
sites. This shows that the (pseudo-) ground state of our 
system is indeed highly degenerate: it is enough to move 
around only a few electrons in order to get into another 
energetically close local minimum. This well illustrates 
the ultrametric structure of the configurational space. 

The A aj 3 distribution is also calculated for smaller sys- 
tems with the same filling factor. For system sizes above 
ss 25 2 the fraction of the electrons different between min- 
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FIG. 3: Distribution of the difference of the local minima 
A a/ 3 (0 for identical states, 1 for complementary ones). The 
distribution is fitted by a Gaussian curve, with mean at 0.054 
and width 0.007. I0 7 overlap integrals are calculated. The 
inset shows the dependence of the mean and width of the 
distribution on system size. 
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FIG. 4: Energy relaxation of the system energy for two differ- 
ent temperatures. In the two upper plots logarithmic fits are 
shown. At lower temperatures (T = 0.01) the initial "flat" 
region is due to a very long equilibration phase of the system 
to reach the steady state from the local minimum. Even after 
averaging over several time series the fluctuations at higher 
temperatures are still large. At even higher temperatures the 
fluctuations dominate and the logarithmic behavior cannot be 
recognized since due to computational limits no more averag- 
ing can be done in reasonable time. The lower graph shows 
the energy relaxation for the case when the Coulomb interac- 
tion range is artificially limited to 3 lattice sites at T = 0.05. 
One sees (I) that the behavior deviates from the log behavior 
and (2) the amplitude of the fluctuations are much smaller as 
for the case with full Coulomb interaction. 

ima and the width of the distribution stay essentially the 
same. Both increase upon lowering the system size; see 
inset in Fig. 03 For very small systems (below 15 2 ) the 
distribution is not Gaussian anymore since there is only 
a small number of possible local minima. This indicates 
that small systems do not have an ultrametric structure 
and therefore do not represent physical systems. 

Energy relaxation. — Next we analyze the relaxation 
of the system at finite temperatures. Due to the ultra- 
metric structure of the system the relaxation can be de- 
scribed as diffusion or random walk on the tree repre- 
senting the structure 0] (Fig. [[]). Starting from a local 
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minimum, nearest neighbor hopping with Boltzmann fac- 
tors e - Ae ij/ T i s applied to all particles in the system - 
random permutations are used to address the lattice po- 
sitions in one time step. The typical energy difference 
Agy of two adjacent lattice sites i and j is of order one 
in the used dimensionless units. 

In Fig. |4] the time evolution of the system energy is 
shown for two different low temperatures. Since the re- 
laxation is very slow and the fluctuations are significant 
(even after averaging over several time series) and fur- 
thermore increase with temperature, a conclusive dis- 
crimination of the behavior is not possible. However, 
it can be well fitted by a logarithmic law in agreement 
with the result of 0. This lends an additional support 
to the hypothesis of the glassy nature of long-range inter- 
acting, disordered Coulomb systems. Previously the time 
evolution in the Coulomb lattice model was studied via 
Monte-Carlo methods, see e.g. [3 Ei, In the 
energy relaxation observed via evolution of the system 
after adding a few electrons showed 1/f-behavior. 

Finally, we investigated the relation between the re- 
laxation of the system and electron correlations reducing 
artificially the screening length of the Coulomb interac- 
tion (in all other simulations the range is of order of the 
system size). In this case the relaxation deviates from the 
logf law (see lower curve in Fig. [4]) indicating the crucial 
role of long-range interactions in forming Coulomb glass. 

Discussion and conclusions. — In conclusion, we have 
derived that in a system endowed with an ultrametric 
configurational space the local energy minima follow ex- 
treme value statistics resulting in an exponential distri- 
bution of the low-lying energy states. Since in the most 
strongly correlated disordered systems the energy land- 
scape is expected to be characterized by only one energy 
parameter, one can assume that barriers separating the 
low-lying states also obey the same statistics. This sug- 
gests that the nonlinear creep dynamic response and the 
related 1// noise are generic properties of all strongly 
correlated disordered systems. 

Our numerical results for a disordered system with 
long-range Coulomb interaction - including the expo- 
nential distribution of the low-lying energy states and 
the logarithmic time relaxation of the system towards its 
ground state - offer strong evidence towards the existence 
of a Coulomb or electronic glass and relates the formation 
of the glassiness with long-range correlations in Coulomb 
systems. This further explains why pronounced glassy 
properties like aging and memory effects as well as appar- 
ently unbounded 1// noise are observed in electronic sys- 
tems with profound long-range electron-electron interac- 
tions. Of course the specific manifestation of the energy 
statistics in the observable quantities requires concrete 
consideration of the peculiarities of the physical system 



involved and goes beyond the scope of this work. In par- 
ticular, one of the most intriguing questions concerning 
the behavior and the nature of the noise in various dis- 
ordered electronic systems will be the subject of a forth- 
coming publication. 
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